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In this paper, we study the homogeneous one-dimensional bosonic gas interacting via a repulsive 
contact potential by using the improved Gaussian approximation. We obtain the gapless excitation 
1 spectrum of Bogoliubov mode. Our result is in good agreement with the exact numerical calculation 

based on the Bethe ansatz. We speculate that the improved Gaussian approximation could be a 
f"^ ■ quantitatively good approximation for higher dimensional systems. 
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I. INTRODUCTION 



Since the concept of Bose-Einstein condensation (BEC) was originally put forward by Bosc and Einstein, the dilute 
Bose gas, as a many-body system which displays macroscopic quantum phenomena such as superfluidity, has been 
extensively studied theoretically. The microscopic description of BEC started with Bogoliubov theory [H41], in which 
(3JQf the destruction and creation operators for the macroscopically-occupicd lowest-energy mode is specially treated as 
c numbers, known as Bogoliubov replacement. Based on Bogoliubov replacement, the Green's function methods 
were applied to a dilute Bose gas at zero temperature 043 ■ Hugenholtz and Pines Q showed that for a repulsive 
interaction, the pole of the one-particle Green's function, approaches zero for zero momentum, which means a gapless 
^ excitation spectrum (usually we call it as Goldstone theorem 0). P.C. Hohenberg and P.C. Martin described BEC as 
spontaneous global U(l) symmetry breaking by introducing external sources, which are set negligibly small in the end 
Q. The interpretation of BEC as symmetry breaking makes the quantum field-theoretic treatment very convenient, 
in which the expectation value of the field operator describes the density as well as the wavefunction of the condensed 
bosons and hence is also called " macroscopic wavefunction" . The effective action approach fiol - flij is usually employed 

■ and kinds of approximations can be easily formulated in this framework, such as Bogoliubov approximation, Pop ov 
£h \ approximation and Hartree-Fock Bogoliubov (HFB) approximation as discussed in detail in the references |15l - l20j . 

However for the bosonic model, if we use the simplest non-perturbative calculation, Hartree-Fock Bogoliubov 
i ^ , ' (HFB) approximation, the spectrum obtained is gapped even in the broken phase. Goldstone theorem is violated 
in such approximation 0, [l5(. Though in the Popov approximation, the spectrum remains gapless, the method is 
t— I ■ not self-consistent and we will also show that we can not apply this method to one dimensional bosonic model. $ 
^ , derivable theory, self-consistent approximation method beyond HFB, including some higher two particle irreducible 
*^ ■ (2PI) diagrams to the effective action, is often used in studying BEC systems. The spectrum obtained in the <3> 
| derivable theory is also gapped [2l| . 

■ In the self-consistent theories such as Hartree-Fock Bogoliubov (HFB) approximation, the Ward identity from 
U(l) symmetry is not preserved due to partial resummations of some Feymann diagrams. Therefore, the Goldstone 
theorem is violated and the resulting excitation spectrum is gapped even in the symmetry breaking phase. In order 
to preserve the Ward identity, we should incorporate the contributions of some other Feymann diagrams and thereby 
remove the gap [2ll423T ] . It is called " covariant Gaussian approximation" in [22[ and we will call it " improved Gaussian 
approximation" (IGA). In principle we can apply similar method to the $ derivable theory beyond HFB (we will call it 
the improved $ derivable theory, or IDT in short), but the theory becomes too complex (involving integral equations 

■ which can not be solved analytically) [2l| . 

$_i | In recent years, interest in ID Bose gas has been revived due to its experimental realization with ultracold bosonic 
atoms (23 - f27j . In one dimension (ID), at finite temperature, the excitation spectra are gapped. However, the ID 
Bosc gas at zero temperature contains gapless spectra and the system is algebraic long range order. In a trapped ID 
gas, the Bose-Einstein condensation (BEC) regimes of a true condensate, quasicondensate regime and the regime of 
a trapped Tonks gas (gas of impenetrable bosons) at finite temperature have been identified in [281 ] . The stability 
and phase coherence of trapped ID Bose gases was studied in [29j . Most of the other relevant works are summarized 
in the review article [30| |. In highly anisotropic traps, where the axial motion of the atoms is weakly confined while 
the radial motion is frozen by the tight transverse confinement, the shape of the Bose-condensed systems reduces 
to one dimension. If the characteristic range of the interatomic potential is much smaller than the typical length of 
the radial extension, the system can be described by the Lieb-Liniger model [H], [32[ , in which the contact potential 
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be exactly solved by the Bethe ansatz and two types of excitations (named Type I and Type II) have been found. 



strength g\ rj is given by gio = — „ t a 1D with a\o being the ID scattering length|33l l34j. The Lieb-Liniger model can 
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Type I excitations are gaplcss with a linear dispersion in the long wavelength limit and reduce to the Bogoliubov 
excitations in the weak coupling limit. Type II excitations, the Fermionic excitations which are prominent in the 
strong coupling regime, have no equivalent in the Bogoliubov theory. J. S. Caux et aZ.[35[ studied the one-particle 
dynamical correlation function of the Lieb-Liniger model by using the ABACUS method [HI, for a wide range of 
values of the interaction parameter. 

In this paper we will apply IGA to the ID Bose gas at zero temperature. This system can be described by the 
Lieb-Liniger model (LLM), which has been exactly solved by the Bethe ansatz. We can compare the result of the IGA 
method with the exact one in order to test the precision and validity of the IGA method. In the future, we shall apply 
IGA to 2D or 3D Bose gas at finite temperature, as in high dimension we can not apply the Bethe ansatz method to 
obtain the exact solution, IGA or IDT is the only approach we can rely on. In higher dimension, the quantum and 
thermal fluctuations are weaker than in ID, the result obtained by IGA or IDT should be better qualitatively and 
quantitatively than that in ID. 

In this paper, we shall study LLM by using IGA, and focus our attention on the excitation spectrum. We will follow 
[22I ] and present IGA method by solving Dyson- Schwinger equations which are generated by functional differentiation 
of the effective action. 

We will show that only the Bogoliubov excitation spectrum (or Type I excitation) can be obtained by IGA. By 
comparing with the results of the Bogoliubov approximation and Type I excitation based on the exact solution (3~i1.l35j|. 
we find that the spectrum obtained in this way is good improvement to the spectrum in the Bogoliubov approximation. 
In order to obtain Type II excitation, we speculate that we shall use more general $ derivable theory beyond IGA 
(we will leave it as our future work) . If we study high dimension Bosonic system, there will be no Type II excitation, 
IGA will give more accurate results quantitatively. 

The rest of the paper is organized as follows. In section II we review the basic formulation of one particle irre- 
ducible (1PI) effective action theory and the Dyson-Schwinger equations. We also present the ID bosonic model 
and the Dyson-Schwinger equations for ID bosonic model in this section. In section III we review the traditional 
approximations, such as Bogoliubov approximation, HFB approximation and Popov approximation. In section IV, 
we present improved Gaussian approximation and obtain an improved gapless excitation spectrum. In section V we 
make a comparison with the exact solution of the ID bosonic model [3ll . [35j . Finally, we give a summary and the 
conclusions. Wc put h = ks = 1 throughout the paper with ks the Boltzmann constant. 



II. THE DYSON-SCHWINGER EQUATIONS FOR ID BOSONIC MODEL 



In this section we shall present the general formulations and the model, and set up all the notations and definitions. 
We shall start with the thermodynamic partition function and set the temperature to zero in the end. For a bosonic 
system, the grand canonical partition function takes the form J37j 



Z = J £>[?/;* Vle-^*^ 1 (1) 

with the classical action S[ip* , ip] given by 



dr J d D x(^*d T ij- /i^*V + H[^*,^]) (2) 

where /3 = , /i is the chemical potential and 7-L is the Hamiltonian density, D is the dimension of position 

space (the formulation is valid for arbitary D, however in this paper, we will only carry out calculations for ID). In 
order to obtain the correlation functions of field operators, a generating functional is defined by coupling fields to an 
external source, 



Z[J*, J] = J V^* ^]e~ (s[r ^ ]+J ^ +jr \ (3) 

where J*t/j is a shorthand for J Q dr j d D xJ*(x 1 t)?/>(x, t) and similarly for J-0*. The connected generating functional 
is defined as 



W[J*,J] = - \nZ[J*,J}. 



(4) 
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The one-point expectation value of the field operators can be obtained by the derivatives of the generating functional 
with respect to the external source, 

8W[J*,J] 
6W[J*,J] 

v (x ' r) = T7(x^r (5) 

where y(x, r) = (tp(x, r)), (p*(x,r) = ("(/>* (x, t)) with 



(■••> = / V, V>] ■ • • e-w>w+ J 'W). (6 ) 



Successive derivatives generate multi-point correlation functions, for instance, 

S 2 W 



6J(x)SJ*(y) 



(P(x)M)) c (7) 



where x = (x,t), y = (y, t') and the connected Green's function (ip* {x)ip(y)) = (ijj* (x)ip(y)) — (ip*(x)} (ip(y)). For 
notation compactness, we define 

(j, j*) = (J x , j 2 ),(v^) = (^i,V 2 ),(^) = W 2 ), 

G mn (x,y) = (ip m (x)ip n (y)) cl (8) 
where m = 1, 2, n = 1,2. G mn (x, y) is related to VF[J*, J] by the following equation, 

G mn (x,y) = . (9) 

dJ m (x)dJ n (y) 

The 1PI effective action is defined by the Legcndrc transformation, 

r[<p*,ip] = w[j*,J}- J?*, (io) 

which is a functional of the field expectation yj*and ip. In analogy with Eq.([5]) the external source can be obtained 
by the derivatives of the effective action with respect to the one-point expectation of the field operators, 

6T[iP%lfi] -J*(x,t), 



<V(x,t) 
ST[<P'M 



-J(x,r). (11) 



V(x,t) 

£i"vr Tim^ri r~\~n o I Icmrr rno rihnin vn In iin l/in Into 



The effective action is the generating functional for vertex functions. Using the chain rule to calculate 7^7^ , we 
have 



6(p m (x) _ y> I" dz S( Pm{x) SJi(z) 



&V>n{y) J SJi(z) 5(p n (y) 



On the other hand, 



8<p n (y) 

Thus by combining Eqs. (|T!?))([T3")) one obtains 



r 5 2 W S 2 T ^ 



S(p m {x) 

-5 mn 5(x-y). (13) 



^ f dzG mi (x,z)T in (z,y) = 5 mn 6(x - y) (14) 
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where T mn (x, y) = s ^ ^g^f 2 ^ and G mn (x, y) is defined in Eq.©. The 1PI effective action T[ip* , tp] can be approxi- 
mately obtained by loop expansion (38l |. 

Dyson-Schwinger equations can be obtained by using the following identity, 

-(s^^+JVW) = 0j (15) 

oip*{x) 



which leads to 



+ J(x) = 0. (16) 



Derivatives of Eq.(fT6|) with respect to the average field tp m {x) shall produce a series of Dyson-Schwinger equations, 
such as 



6 / 5S[P,if>] \ | S 
8<p(y) \ Sip*(x) J 5(p(y) 



5(p[y) 

Successive functional derivatives with respect to tp (z) yield higher order Dyson-Schwinger equations, which involve the 
correlation functions of more field operators. Therefore, the infinite Dyson-Schwinger equations must be truncated to 
form a set of closed equations in order to carry out any calculations. Let us term Eq. (|16|) as the first Dyson-Schwinger 
equation and Eq. (|17|) as the second Dyson-Schwinger equation. 

We apply the Dyson-Schwinger formalism to a system of one-dimensional bosonic gas interacting via a repulsive 
contact potential, described by the Lieb-Linigcr Hamiltonian 

N N 

H = -^(df i /dxl)+gJ^S{x i -x j ), (18) 

i— 1 i<j 

where the mass of the particle has been set to 2m = 1 and g is the contact interaction strength, which is related to 
the ID scattering length experimentally. The second quantization form reads 

H= Jd D X (V (x)(-V 2 Wx) + ^ to^M^to^to) , (19) 

where we have used the notation for a general position space dimension D and bear in mind that we will study the 
ID case of D = 1 in the end. 

In path-integral formalism, the grand canonical partition function takes the form 

Z = J VU>*i>}er s ^*^ (20) 

with the classical action S[tp* , ip] given by 

dr J d D x (V {d T -H- V 2 ) tp + ^-gip*ip*if}i?J (21) 
where ip = ip(x, r), f3 — -r-Kp and [i is the chemical potential. By variable rescaling 

x = 5 - 1 x', M = . 9 2 /i', (22) 
the action can be recast as a simple form dependent only on one parameter \J , 

^ dr' J d D x! (V* (dr- - V 2 , - //) 1/ + ^'VW) • (23) 
In the following discussions, we will omit the primes for simplicity, 

S [ip*,ip] = ' dr J d D x (V {d T -S7 2 -pi)i>+ Jff#J . (24) 
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Starting with the rescaled action in Eq. ([24|) , we define the generating functional 



Z[J*,J] = J Viij^^e-^'^+J^+^'l 
The first Dyson-Schwinger equations take the form 



(25) 



(d T - v 2 - n) <p 2 + (ihihih) + Ji = o, 

(-d T - V 2 - fi) <p X + (Vi^i^) + J 2 = 0, 
where implicitly all the arguments are x = (x, r) . By Wick theorem we know 

where (•••)„ means connected correlation functions. Substituting Eq. (f2~Tj) into Eq. 



yields 



(26) 
(27) 



(d T - V 2 - ft) (f2 + ifiipl + <^iG 22 + 2ip 2 Gi2 + (ipifafa) c + J\ =0, 
(-<9 r - V 2 - /i) <pi + iplip 2 + tp^Gu + 2^iGi2 + {ipzipiipi) c + J 2 = 0, 



(28) 



where all the default arguments are x = (x, r) and Gu = Gn(x,x) , G22 = ^22(^5^), G12 = (^i{ x )^2{%)) c - Gij — 
Gij\X)X} is a constant for a translational symmetric system which is the case in this paper. Further differentiations 
of Eq. (|2"5]) with respect to fi{y) and <f2(y) result in the second Dyson-Schwinger equations, 

rn(x, y) = ((fij + G22) 5(x - y) 

+^iA 22 i(a;, y) + 2(p 2 A 12 i{x, y) + - — — (ijjifofo) c , 



8<pi(y) 

^22(x,y) = (ipf + Gn) S(x - y) 



+^ 2 Aii 2 (a;, y) + 2ipiA 12 2{x, y) + 



<Wy) 



(ip 2 ipiipi) c , 



Ti 2 (x, y) = (d T - V 2 - fi + 2ip x y 2 + 2G 12 ) S(x - y) 



+ipiA 222 (a;, y) + 2tp 2 A 122 (x, y) + 



(^1^2 ^2)0 



T21 (x, y) = (-d r - V 2 - fi + 2^2 + 2Gi 2 ) S(x - y) 
+(p 2 Am(x, y) + 2cp 1 Ai 2 i(x, y) + - — i— {ifaipiipi) c , 



<Vi(y) 



(29) 



where x = (x, r), y = (y, r') and A m „;(a;, y) = l5G ^"^' x - ) with m (n, Z) = 1, 2. Since what we consider is a homogeneous 
gas, we can set 



<pi(x,r) = V32(x,t) = u, 
where u is a real constant number. Further, we define the Fourier transformations 



(30) 



5(x — y) = 
A mn i(x,y) = 



duj r d D k 

2^ J (2tt) d 
d D k 



G mn (x,y) — I — 



(2tt) d 
d D k 

D 



2tt J (2?r) 



d D k 



2tt J (2tt 



ik-(x— y)-io)(r— t') 

A mni (fc)e ik - (x - y) -^ (T - T,) , 
G„ m (fc)e lk - (x - y) -' i " (T - T ' ) , 
r m P '' k '( x ~ y )~'"( T ~ T ') 



(31) 



where k = (k, cj) and w denotes the Matsubara frequency in the zero temperature limit. In the frequency space, 
Eq. dTil) is recast as 
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G im (k)T mj (k)=S ij . (32) 

m=l,2 

The first and second Dyson-Schwingcr equations arc not closed equations. They are impossible to solve unless 
truncations are performed. 

III. THE TRADITIONAL APPROXIMATIONS 

The traditional approximations, such as Bogoliubov approximation, HFB approximation and Popov approximation, 
have been exhaustively discussed in the literature. In order to clarify the interrelations of the various familiar schemes 
and the IGA scheme we shall present later, in this section we formulate those approximations by truncating the first 
and second Dyson-Schwingcr equations. 

a. Bogoliubov approximation: Ignoring any correlations, only the first Dyson- Schwinger equations Eq. (|28[) are 
retained: 

(d T - V 2 - n) v? 2 + fifl + J i = °> 
(-Or - V 2 - fl) tpi + ^2 + h = 0, (33) 

and the two-point vertex functions are defined by Tij{x, y) = — fy'.^j \ j i { x )=o where Ji(x), fj(y) are related by Eq. ([33| . 

rii(x,j/) = Lp\5{x - y), 

T 2 2(x,y) = (pfS(x-y), 

Ti2(x,y) = (d T - V 2 - n + 2Lpiip 2 ) S(x - y), 

T 21 (x,y) = (-5 T -V 2 -/i + 2<pnp3)S(x-y). (34) 

By using the homogeneous and static condition in Eq. (|30[) and applying the Fourier transformation in Eg. pip , wc 
rewrite Eq. (|3"3")) when Ji(x) = as 



v z = n (35) 



and Eq.(|34]) becomes when Ji(x) = 0, 



r n (fc) = v\v 22 {k) = v 2 , 

Ti 2 (fc) = -iw + k 2 +w 2 , 

r 2 i(fc) = ico + k 2 +v 2 . (36) 
With the help of Eq. (j3"2")) we obtain the Green's functions in Bogoliubov approximation, 

G u (fc) Gia(fc) \ = 1 

G 2 i(k) G 22 (fc) J (i w )2 _ k2 (k 2 + 2v 2 ) 

--(JJ (37) 
The Bogoliubov spectrum is given by the pole of the determinant of Matrix Eq. (|37[) 



e B o g (fc) = fc v / fc 2 + 2u 2 . (38) 
In this approximation, the particle density n is equal to v 2 . 
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b. HFB approximation: If two-point correlation functions are kept, ignoring three or higher point correlation 
functions, the first Dyson-Schwinger equations Eq. (j2"5)) become 

[pr - V 2 - m) ¥>2 + + V1G22 + 2(p 2 G 12 = 0, 
(-dr - V 2 - (i) <f! + w\<p 2 + f2G n + 2<^iGi2 = 0, (39) 



and the second Dyson-Schwinger equations Eq. ([29|) become 

Fii(x,y) = (tp\ + G 22 ) S(x - y), 

F 22 {x,y) = (tpi + Gn) S(x - y), 

T 12 (x,y) = (<9 T -V 2 - f i + 2ip l (p 2 + 2G 12 )5(x-y), 

r 2 i{x,y) = (-d T -Vl-^ + 2^2 + 2G 12 )5(x~y). (40) 

By using the homogeneous and static condition in Eq. (|30l) and applying the Fourier transformation in Eg. pip , we 
rewrite Eq. (|3T))) as 



= -/x + v 2 + Gn + 2Gi 2 , Gn = G 22 , (41) 



and Eg. ([iOll as 



ru (A) - v 2 + G n , 

r 22 (k) = v 2 + G lu 

T 12 (k) = -icu + k 2 +v 2 -G u , 

r 2i (fc) = iuj + k 2 +v 2 -G n - (42) 
With the help of Eq. (|32|) we obtain the two-point Green's functions in HFB approximation, 

f Gn(fc) G 12 (k) \ 1 

^ G 2i (fc) G 22 (fc) J (tujy - (k 2 + 2v 2 ) (k 2 - 2G n ) 

/ v 2 + Gn iu - (k 2 + v 2 - Gn 

X V -*w - (k 2 + v 2 - Gn) v 2 + G n 

Then the HFB spectrum is given by 



(43) 



£HF B (k) = V(k 2 +2u 2 )(k 2 -2Gn). (44) 
The variable Gn can be determined in a self-consistent way. By the definitions of Gnand Gi 2 , there are 

1 f°° 1 

Gn = --r- (v 2 + Gn) / dk - (45) 

4^ v 1 J_ x x /(fc 2 + 2u 2 )(fc 2 -2Gn)' V ' 

and 

1 f 00 ( (fc 2 + u 2 -Gn) \ , x 

Gi 2 = — / dk - V 11 J - 1 . 46 

4^7-oo VV(fc 2 + 2w 2 )(fc 2 -2Gn) ; 

In HFB approximation, the particle number density is 

n = v 2 + Gi 2 (47) 
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c. Popov approximation: Popov approximation is well-known for its gaplcss excitation spectrum. It differs 
from the HFB approximation in neglecting the "anomalous" two-point correlations Gnand G22, so that the Dyson- 
Schwinger equations take the form 

-H + v 2 + 2G 12 = (48) 

and 

r u (fc) = v 2 ,T 22 (k) = v 2 , 

T 12 (k) = -iiv + k 2 +v 2 ,T 21 (k) =iLo + k 2 +v 2 . (49) 

In terms of the variable u 2 , the two-point Green's functions have the similar form as those in the Bogoliubov approx- 
imation, 

Gn(fc) G 12 (k) \ 1 



G 21 {k) G 22 (k) J («j)2 _ k 2 ( k 2 + 2u 2) 

(50) 



e + v 2 ) v 2 



and also the excitation spectrum 



By the definition of G12, there is 



epopov(fc) = *Vfc 2 + 2iA (51) 



fdu f d D k l u J -(k 2 + v 2 ) 
12 J 2ir J {2ir) D (iw)2 - k 2 (k 2 + 2v 2 ) ' 1 j 

The particle number density is given by n — v 2 + G\ 2 . However, in ID, the above equation leads to 

n = v 2 + -1- I" dk (-1 + ^k 2 + 2v 2 . V o = ) , (53) 

V k kVk 2 + 2v 2 J 

which is infrared divergent. So the Popov approximation is inapplicable here. 

The reason for Popov theory to break down in ID is that phase fluctuations are not considered properly. Ref.(39j 
gave a detailed discussion of this problem and proposed the modified Popov theory, in which the inappropriately 
incorporated phase fluctuations are subtracted and thus the infrared divergence is removed. The particle number 
density from the modified Popov theory shall be given by 



00 



n = v2 + T" / * ( -1 + ^fW=^ ) . ( 54 ) 



4W-OC V Vk 2 +2vy ' 

which is free of divergences. 



IV. IMPROVED GAUSSIAN APPROXIMATION 



In this section we shall present another strategy, IGA (improved Gaussian approximation) which takes account of 
quantum fluctuations more precisely (adding some Fcymann diagrams to preserve symmetry requirement) and retains 
the gapless Goldstone mode. 

By preserving up to two-point correlation functions in the first Dyson-Schwinger equations, however we will keep 
source terms here for a while in order to define the Green's function in IGA scheme. 

(d T - V 2 - n) ip 2 + ipiipl + ipiG 2 r 2 + 2(p 2 G l { 2 + Ji = 0, 
(-d T -V 2 - f i)i Pl +^ 2 1 i P2 +^ 2 G t { 1 +2i Pl G t 1 r 2 + J 2 = 0, (55) 

and 
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FIG. 1: Feynman diagrams for the corrections to the two-point vertex function obtained by HFB approximation. 



T\\{x,y) = (tf% + G&)6(x-y) 

T%(x,y) = (a r -Vl-n + 2<p 1 <p 2 + 2G&)6(x-y) 

ffifrv) = {-d T -V* - n + 2<p 1 tp 2 + 2G&)5{x-y), (56) 



where tr is the abbreviation of " truncation" . We will define 

Tij (x, y) = — 5 x Jt< f\ I J»(x)=o (57) 
where the relations between Jj (x) and ipj(y) are given by Eqs. (|55l56l) . 

r u (x,y) = (^ + G^ 2 )(5(x-y) + ^ 1 A t 2 r 21 (a;,y) + 2^ 2 Af 21 (2;,y), 

r 22 (x ;2/ ) = (p? + G£) <5(x-y) + <^ 2 Af 12 (a;,y) + 2^ 1 Af 22 (a;,y), 

r 12 (x, y) = (d T -Vl-fi + 2^ 2 + 2Gf 2 ) S(x-y) + ViA 222 (x, y) + 2ip 2 A? 22 {x, y), 

T 21 (x,y) = (-d T - VI - f i + 2i Pl <p 2 +2G t 1 r 2 )6(x~y) 

+<p 2 A t { 11 (x,y) + 2<p 1 A t { 21 (x,y), (58) 

where 

^n^,y)^ 6G f;^\ (59) 

and in the end we shall take Ji(x) = 0. From Tij(x,y), we can obtain the Green's function which is the inverse of 
Tij(x,y). The result obtained is gapless [H, H(|. T*j(x,y) is obtained from the truncated Dyson-Schwinger equation 
ignoring three-point Green's function. We comment that diagrammatically the corrections ATij(x, y) = Tij(x,y) — 
TYj{x,y) correspond to some additional diagrams (2l| - [23j . which are plotted schematically in Fig. [TJ In the Feynman 
rules of Fig. [TJ the point vertices are defined by the interaction part of S [ipi + ip\ 1 ip 2 + ip 2 ], i.e., the part with three 
and four ipi fields expanded around ifi{. The lines in Fig. [T] stand for the truncated Green's function G*J. The cross 
in Fig. QJepresents (pi (details can be found in Ref . [2l| ) . By using the homogeneous and static condition in Eg. ([50)1 
and applying the Fourier transformation in Eg. pip , we rewrite Eq. (|58[) as 



r n (fc) = v 2 + G t 1 r 1 +vA 2 r 21 (k) + 2vA\ r 21 (k), 

T 22 (fc) = v 2 + G t 1 \+vAT 12 (k) + 2vA^ 2 (k), 

T 12 (fc) = -tuj + k 2 +v 2 -Gf 1 +vA t 2 r 22 (k) + 2vA t 1 r 22 {k), 

r 21 (fc) = lui + k 2 +v 2 -G\ r 1 +vA t 1 r 11 (k) + 2vA t 1 r 21 {k). (60) 

v 2 and G*J in the above equation are obtained from the HFB equations in the previous section as in the end we take 
Ji(x)=0. 

We start to calculate A^ nl (x,y). First, we differentiate Eq. (|5"6"|) with respect to (p m {z), 
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riii(x,2/,z 

(x,y,z 



221 



121 



(x,y,z 

^211 (-^i Vl Z 
-^222 i X i Vi z 



^11 2 {.X, 2/> z 
Ti22 { X l Vi z 
^212 (^j 2/j z 

where z = (z,r"). From Eq.(|14j) we know 



<5(.T 
<5(.T 

8(x 
8(x 
8(x 
8(x 
8(x 



y) (A? n (x,x,z) + 2ipi8(x-z)) , 
y) {2A t { 2l {x, x, z) + 2^ 2 8{x - z)) , 
y) (2A? 21 {x,x,z) + 2<p 2 6(x-z)), 
y)&n 2 (x,x,z), 

y) ( A 222 fal Z ) + 1ip 2 8(x - Z)) , 

y) (2A\ r 22 (x, x, z) + 2ip 1 8(x - z)) , 
y) (2Af 22 (x, x, z) + 2<p 1 8(x - z)) , 



(61) 



(62) 



The derivatives of the above equation with respect to ipi (z) result in 



a tr 

rani 



(x,y,z) 



Y.T.j dx ' j dy'Gt m ,(x,x')T^ n>l (x',y',z)Gi r , n (y',y). 



(63) 



One can now take Ji{x) = 0. G^ m /(x, x') is thus given by HFB approximation in the above equation. By substituting 
Eq. ([6T|) into Eq.([63|) and setting x = y, one obtains a set of closed equations for A mn ;(x, x, z), 



A mnl {x,x,z) = -2v z)G£(z, a:) + G£^ 

-2 J dx' (G*,(x, x')Gt(x\ x)A m (x', x', z) + G^x, x')Gt(x', x)A m (x\ x', z)) 

- j dx'(G t r i (x,x')G*( x >, x )A m (x\x\z) + G t rfa (64) 

where I is defined by 8^ = 0, which means / = 1, I = 2 or I = 2, I = 1. By applying the Fourier transformations in 
Eq.(j3l|) we rewrite Eq. ([64]) as 



A^(fc) = AT u {k)I m u n {k) + A%{k)I ml , ln {k) + A%{k) (2I mlJn (k) + 2I mlln (k)) 

+2^ (l mlJn (k) + I mTtln (k) + I mIJn (k)) , (65) 

where 

Imn,m'n'(k) = - j ^ J ^G^(fc X + k)G^ n , (h) (66) 

and the two-point functions G^ n (fc) are those obtained from the HFB equations. We can explicitly integrate u>\ in 
Eq.(|551). for example, 



1 f d D ^ -(^ 2 + Gf 1 ) 2 



n ' U{ ij {2n) D ^/p + ki) 2 + 2u 2 ) ((k + kO 2 - 2G\\) (k 2 + 2u 2 )(k 2 - 2G\\) 
x( - 

\uj + y/((k + kO 2 + 2v 2 ) ((k + kj 2 - 2Gf 1 ) + ^(k 2 + 2u 2 )(k 2 - 2G\\) 

^/[{k + k^ + 2v 2 ) {(k + k x ) 2 - 2G\\) - ^(k 2 + 2u 2 )(k 2 - 2Gf x ) ^ (67) 
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which shall be used for analytic continuation described below. Next, we insert the A mn i(k) solved from Eq. (|65j) into 
Eq.([60p. so as to obtain the improved two-point vertices T mn (k). The improved two-point correlation functions take 
the form 

G 22 (fc) = [v 2 + G\\ + v (A*/ 21 (fc) + 2A*' 21 (fc))] , 

Gn(fc) = [v 2 + G\\ + v (K\\ 2 (k) + 2Af 22 (fc))] , 

G2l(fc) = w(kj ^ lbJ (k2 +v "~ G " + v (Arn(fc) + 2Ar2l(fc)))] ' (68) 

where M(k) is the determinant of the matrix T mn (k) and reads 

M(k) = [v 2 + G% + v (Ar 21 (fc) + 2Af 21 (£;))] 
x [v 2 + G\\ + v (Af 12 (fc) + 2Af 22 (fc))] 
- [-uj - » + k 2 + 2v 2 + 2Gf 2 + v (A 2 r 22 (k) + 2A\ r 22 (k))] 

x [%u - ft + k 2 + 2v 2 + 2G\ r 2 + v (A? n {k) + 2Af 21 (fc))] . (69) 

The Green function in Eq. (|68p gives a gapless excitation spectrum, which shall be shown by the numerical result 
and also can be analytically verified by investigating the poles of the Green's function. Analytically, one can prove 
M(0) = to make sure that the excitation spectrum is gapless. The details of the proof are put in the Appendix. 

One can obtain the real time Green's function, retarded and advanced Green's function by analytic continuation, 
iuj — > fl ± irj, where r\ is an infinitesimal positive number. The spectral weight function is then obtained by using the 
relation [4l[ 

p (k, ft) = 2lmG A (k, fi) = -2ImG* (k, 0) . (70) 

Eq. (|67[) is an analytic function of "complex" variable except on the real axis in f2 plane. The retarded and 
advanced Green's function obtained therefore have desirable analytic properties. 

There is an equivalent formalism of the IGA approximation in the framework of the improved <f> derivable theory. 
The $ derivable theory can start with the two particle irreducible (2PI) action functional V [ipi,ip 2 ,G] which takes 
the form 



T[ip u cp 2 ,G] = S [cp u cp 2 ] + \Tr In G 1 + ^Tr [B' 1 (G - D)] +9[<p 1 ,<p a ,G] (71) 



D 



-i\ S 2 S[(pi,ip 2 ] 



where ((pi,ip 2 ) = (¥>*,¥>) as defined previously, and G represents matrix (G)^ = Gij of Green's functions. In the 
order of HFB approximation (omitting higher order diagrams like the setting sun diagram), 

$ [(P! t <f2, G} = ^Jdx [Gn (x, x) G 22 (x, x) + 2G\ 2 (x, x) G 2 i (x, x)] (72) 

We will obtain the same equations as Eq. (|39|) and Eq. (|40| of the HFB approximation if we require 



ST [yi,y? 2 ,G] = Q ST [y?i,y 2 ,G] = Q 
5<fi SGij 

In the framework of the $ derivable theory, IGA can be reformulated as below Ref.[2l[. The IPI effective action 
T[tpi,tp 2 } is equal to f [<pi, ip 2 , G tr (<pi, <p 2 )} with G tr (ipi, (p 2 ) defined by ^^ffi 2 ^ |G=G"-(yi,y 2 ) = 0. Then from 
T[<pi,(p 2 ], one obtains the inverse Green's function Ty = - g^g^f 2 ^ = ^^'g^'.'gy, ^ tpi ' tp2 ^ m F or technical details, 

see Ref. 21|. Substituting the solution of Eq. ([73|) to the functional T [ipi, <p 2 , G], we obtain a quantity T. The 
thermodynamical potential is f3~ 1 T. According to the thermodynamical relation, the particle number density n is 
equal to — with L being the size of the ID system( L is infinity in the thermodynamic limit). Using Eq. (|73|) and 
Eq. ([7Tj) . we know the density n is equal to v 2 + G* 2 , the same as the case of HFB. It is also valid in any $ derivable 
theory or improved $ derivable theory beyond HFB. 
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COMPARISON WITH THE EXACT SOLUTION 



The references 35 [present an exact solution of the Lieb-Liniger model, which gives the exact excitation spectrum. 

The references [31 , |35| consider a one-dimensional system of length L (satisfying periodic boundary conditions) , 
with TV bosonic particles interacting via a repulsive contact potential of strength 2c, governed by the Lieb-Linigcr 
Hamiltonian 



N 

H = -J2(d 2 /dx 2 l )+2c ]T 5{ Xi - Xj ). (74) 

i=l l<i<j<N 

The excitation spectrum is plotted as u /n 2 ~ k/n , with n being the particle number density ~ . The dimensionless 
parameter of the system is defined by 7 = Comparing Eq. ([T8]) with Eq. ([74| , there is g = 2c and hence the 
corresponding parameter in the field-theoretic treatment takes the form 7 = ^.Comparing Eqs. (|21[)(|23[) . we know 
n = gn', k = gk' and w = g 2 uj', which implies uj/n 2 = ui'/n' 2 , k/n = k'/n', where we restore the notation k', n', 
ui' for the rescaled quantities after Eq. ([2"3"|) (we had dropped prime for simplicity). Therefore, in order to compare 
with the exact solution, we should plot the excitation spectrum in the form ui 1 /n' 2 ~ fc'/n' with n' being the rescaled 
particle number density, at the parameter 7 — ^7. At the parameters 7 = 1, 7 = 32, 7 = 64, corresponding to 
n' = 1/2, n' = 1/64, n' = 1/128, we plot the spectrum obtained from the different approximation schemes in Fig. 
[3J The IGA spectrum, which incorporates extra corrections based on the HFB spectrum, is gapless, while the HFB 
spectrum is gapped. When the particle density is high (7 is small) , all approximation schemes lead to good results, 
which implies that quantum fluctuations are weak at a high particle density. Furthermore, at a very low particle 
density when quantum fluctuations become strong, the IGA scheme shows its advantage. Specifically, at 7 = 32 and 
7 = 64, the IGA spectrum is in good agreement with the exact one, while the Bogoliubov spectrum is not accurate 
quantitatively. 



VI. SUMMARY 



We have presented IGA (improved Gaussian approximation) to treat one dimensional bosonic gas. The Green's 
function obtained by IGA satisfies Ward identities from U{\) symmetry and therefore the spectrum is gapless. 

We have formulated all the traditional approximations (Bogoliubov approximation, HFB approximation and Popov 
approximation) in terms of truncations of Dyson-Schwingcr equations. The HFB approximation is the well-known 
self-consistent approximation, but it leads to a gapped excitation spectrum, violating the Goldstone theorem. The 
spectrum obtained by IGA scheme, which incorporates more quantum corrections to the HFB spectrum, is gapless. In 
order to test the validity and precision of the IGA method, we apply it to the one-dimensional bosonic gas described 
by the Lieb-Liniger model. We can only obtain Type I excitation ( Bogoliubov spectrum) within IGA method. 
Nevertheless, by comparison with the type I spectrum exactly solved by the Bethe ansatz, we find that the IGA 
method gives quantitatively good results on type I spectrum. 

The idea of the IGA method can be applied to improve higher order $ derivable theory (the HFB theory is the result 
of the lowest order $ derivable approximation) [2l], [22| ■ The essence of the idea is to add extra Feynman diagrams 
to preserve the symmetry of all the Feynman diagrams, and thereby restore the Ward identity. The IGA method 
makes improvement based on the HFB approximation. When quantum fluctuations are very strong, higher order $ 
derivable approximation beyond the HFB approximation will be required and then the corresponding improvement 
to restore the Ward identity can be performed in a similar way. In order to get type II (Fermionic excitation), one 
probably shall go beyond IGA and use "improved" high order $ derivable theory. 

The IGA method presented here can be employed to handle many other Bose-condensed systems, including 2D 
or 3D, zero temperature or finite temperature, homogeneous or in optical lattices. For higher dimension systems, as 
there is no type II excitation and quantum fluctuations are weaker, IGA shall be expected to give more quantitatively 
accurate results. 

As one of applications of IGA, we have carried out the IGA calculation on type II superconductor where acoustic 
and optical spectra are obtained non-perturbatively. The results will be presented elsewhere [42j . 
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FIG. 2: (color online)The spectra obtained from Bogoliubov approximation (green lines), HFB (blue dashed lines), IGA (red 
dotted lines) and exact numerical calculations (black squares) from Bethe ansatz [35l | are compared at three different 7. 
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VII. APPENDIX 



We shall prove M(0) = 0, namely, 

= [v 2 +G%+v{At 21 (0) + 2Af 21 (0))] 



x [v 2 + G\\ + v (Af 12 (0) + 2Af 22 (0))] 

- [-» + 2v 2 + 2Gf 2 + v (A*£ 2 (0) + 2Af 22 (0))] 

x [-/x + 2t, 2 + 2Gf 2 + v (Af n (0) + 2Af 21 (0))] . (75) 



By using the homogeneous and static condition in Eq. ([30| and applying the Fourier transformation in Eq. (f3T|) . we 
rewrite Eq. ((55|) as 



and Eq. ([55|) as 



= —/j, + v 2 + G*i + 2G* 2 , 
G\\ = (76) 



rfi(fc) = v 2 + G&, 

rr 2 (fc) = v 2 + g{\, 

rf 2 (fc) = -ic + i^ + u 8 -^, 

T&(k) = iw + k a + u a -G&. (77) 



Note that the value of the external sources Ji and J2 has been set to zero. With the help of Eq. (|32| we obtain the 
two-point truncated Green's function, 



G\\{k) Gf 2 (k) \ _ 1 

Gti(k) Gr 2 (fc) J {iu) 2 - (k 2 + 2^ 2 ) (k 2 - 2G\\) 



GW iuj - (k A + v z - Gf x 



(78) 



■ioj - (k 2 + v 2 - v 2 + G\ r 1 

which has the same form as the Green's function in the HFB approximation. Using Eq.(|76|) we can rewrite Eq. (|75|) as 

= [ w 2 +Gr 1+w (Ar 21 (o)+2Ar 21 (o))] 

x [u 2 + Gr i +u(A t 1 'i 2 (0) + 2A t / 22 (0))] 
- [v 2 - G\\ + v (Ar 22 (0) + 2Af 22 (0))] 

x [v 2 - Gf, + v (Af n (0) + 2Af 21 (0))] . (79) 
By inserting Eg. ([75]) in Eg.(|Ml). it is easy to verify that 

/ mM „(0) = / Mife (0), (80) 

where m,n,l,fh,fi,l = 1,2 with the constraint S rn m = 0, 5 n n = and % = 0. For example, 111,11(0) = /22.22(G), 
i"2i,n(0) = /ia,3a(0), etc. Eq.® and Eq.gSJ lead to 



A^(0) = A*T fir (0). (81) 



There is 



A% nl (k) = A t r [ ml (k), (82) 

which is evident from the definition in Eq. (|59")) . Using Eq. (j5Tj) and Eq. (j52")l we can rewrite Eq.(fT9"|) as 

= (2v + 4Af 21 (0) + Af n (0) + A&i(0)) 

x {2G\\ + vA% 21 (0) - vA? n (0)) . (83) 
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From Eq. (|65]) wc find that 



(\ tr (0\ \ tr - 2V ( J H.ll(°) ~ J 12.2l(°)) (OA-. 

(A 221 (0) - A m (0)) - (1 + J l(0 )_j 12i2l( o)) • ( 84 ) 



By straightforward calculations, wc know 

wo)-wo) = ^r.ft 1 . <»> 



and 



1 r°° 1 

Gf = - (v 2 + Gt) — \ dk - (86) 

Eq.® follows from the definition of Gft, that is Gft = Gftfox) = / / J^Gf^k, w). By comparing Eq.([55l) 
and Eq.(|86]) we know 



/u,ii(0)-I 12 ,2i(0)= ,„, 2 , " tr , - (87) 



Eq.JMJ) and Eq. (|87jl lead to 

2Gf 1 + t;A£ 31 (0)-t;AS r u (0) = 0. 
Thus Eq. ([8"3"|) is proved and also Eq. ([75|) is proved. 
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